!------------------------------------------------------------
! $Id: x2y2model.shar,v 1.2 2011/05/14 21:34:35 daw Exp $
! Description: does MD by velocity verlet
! Author:  Murray S. Daw <daw@clemson.edu>
! Created: 6 Nov 2007
!-------------------------------------------------------------
! Philip Ladd contributed some lines in 2006
! as a mod to original by msd (using RK method)
! Velocity verlet is sympletic and time-reversible
!	q(n+1) = q(n) + p(n)dt + 0.5f(q(n))dt^2
!	p(n+1) = p(n) + 0.5(f(q(n)) + f(q(n+1))dt

module md_vv

! Load modules
  use prec
  IMPLICIT NONE

  real(r8) :: lambda               ! lambda (anharmonicity)
  real(r8), dimension(2) :: fs      ! forces
  real(r8), dimension(2) :: fsnext  ! forces (next step)
  real(r8) :: dt, mr ! timestep (set in setupmd1dvv)

  private :: fs, fsnext, dt, lambda

contains

  subroutine setupmdvv( lambdain, dtin )

    real(r8), intent(in) :: lambdain
    real(r8), intent(in) :: dtin

! store lambda
    lambda = lambdain
    dt = dtin

    return
  end subroutine setupmdvv

! ps is momentum NOT velocity
  subroutine mdvv( qs, ps, nsteps )

    real(r8), dimension(2), intent(inout) :: qs
    real(r8), dimension(2), intent(inout) :: ps
!    real(r8), intent(in):: mass
    integer, intent(in) :: nsteps
    real(r8):: radius, midvalue, e0, ef
    real(r8), dimension(2) ::qsn, psn, pupdate, bounce


! Declare local parameters
    integer :: istep        ! index for summations

!   mr = mass
    radius = 1.
    bounce = 0.
! calculate forces in cartesian coordinate

    fs(1) = -(qs(1)-2*lambda*qs(1)*qs(2))

    fs(2) = -(qs(2)-lambda*(qs(1)**2-qs(2)**2))
    
!	e0 = (qs(1)**2 + qs(2)**2 + ps(1)**2 + ps(2)**2)/2. - lambda*(qs(2)*qs(1)**2-qs(2)**3/3.)
! Numerically integrate using Velocity Verlet
    do istep = 1,nsteps


       qsn(1) = qs(1) + ps(1) * dt + 0.5*fs(1) * dt**2
       qsn(2) = qs(2) + (ps(2)) * dt + 0.5*(fs(2)) * dt**2

! H-H potential only applies within the circle qs(1)**2+qs(2)**2=radius
! particle is also bounded within the same circle

	If ((qsn(1)**2+qsn(2)**2) .gt. radius**2) then

		bounce(1) = -(dt*qs(2)*(qs(1)*ps(2)-qs(2)*ps(1)))/(qs(1)**2+qs(2)**2)
		bounce(2) =  (dt*qs(1)*(qs(1)*ps(2)-qs(2)*ps(1)))/(qs(1)**2+qs(2)**2)
			If(bounce(1)*ps(1)+bounce(2)*ps(2) .lt. 0) then
				bounce = -bounce
			end if			
		qsn = qs + bounce


!		psn(1) = (-qs(1)**2*ps(1)+qs(2)**2*ps(1)-2*qs(1)*qs(2)*ps(2))/(qs(1)**2+qs(2)**2)
!		psn(2) = (-2*qs(1)*qs(2)*ps(1)+qs(1)**2*ps(2)-qs(2)**2*ps(2))/(qs(1)**2+qs(2)**2)

		psn(1) = (-qs(1)*ps(1)*qsn(1)+qs(2)*ps(1)*qsn(2)-qs(1)*qsn(2)*ps(2)-qs(1)*qsn(2)*ps(2))/(qs(1)**2+qs(2)**2)
		psn(2) = (-qs(1)*qsn(2)*ps(1)-qsn(1)*qs(2)*ps(1)+qs(1)*ps(2)*qsn(1)-qs(2)*ps(2)*qsn(2))/(qs(1)**2+qs(2)**2)
	
		print *, " bouncing, R^2= ", qsn(1)**2+qsn(2)**2
	!	print *, " bouncing, change of kinetic energy: ", psn(1)**2+psn(2)**2-ps(1)**2-ps(2)**2
		ps=psn
	
		qs = qsn
		
                fsnext(1) = -(qs(1)-2*lambda*qs(1)*qs(2))
                fsnext(2) = -(qs(2)-lambda*(qs(1)**2-qs(2)**2))
		fs = fsnext
		print *, " ps(1) = ", ps(1)
	else
	
		qs = qsn
  
        	fsnext(1) = -(qs(1)-2*lambda*qs(1)*qs(2))
		fsnext(2) = -(qs(2)-lambda*(qs(1)**2-qs(2)**2))
             
		pupdate(1) = 0.5*(fs(1)+fsnext(1))
		pupdate(2) = 0.5*(fs(2)+fsnext(2))
		ps(1) = ps(1) + pupdate(1)*dt
		ps(2) = ps(2) + pupdate(2)*dt
		fs = fsnext
	end if
    enddo
!      ef = (qs(1)**2 + qs(2)**2 + ps(1)**2 + ps(2)**2)/2. - lambda*(qs(2)*qs(1)**2-qs(2)**3/3.)

!	check energy conservation
!	print *, "change of energy during MD:", ef-e0
    return
  end subroutine mdvv

end module md_vv
